Thermodynamics of an interacting trapped Bose-Einstein gas 
in the classical field approximation 
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We present a convenient technique describing the condensate in dynamical equilibrium with the thermal 
cloud, at temperatures close to the critical one. We show that the whole isolated system may be viewed as 
a single classical field undergoing nonlinear dynamics leading to a steady state. In our procedure it is the 
observation process and the finite detection time that allow for splitting the system into the condensate and the 
thermal cloud. 
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Successful experimental realization of Bose-Einstein con- 
densation in a dilute gas of alkali atoms [Jj]] offers a new and 
fascinating tool to probe the border between quantum and 
classical worlds. The theory of the dynamical behavior of 
such a many body system is very difficult and cannot be solved 
exactly. It is particularly hard to extract quantitative predic- 
tions in the vicinity of the critical temperature. At the other 
extreme, at zero temperature, the weakly interacting Bose gas 
is well described by the mean-field approach. All particles oc- 
cupy the same quantum state whose wave function is the low- 
est energy solution of the Gross-Pi taevskii (GP) equation [^]. 
At low temperatures the Bogoliubov approximation comes 
handy ^ with its quasi-particles that have just been observed 
in a direct experiment [Eh. The theory gets much more com- 
plex at higher temperatures. There is a considerable effort to 
develop a working theoretical and numerical tool valid there. 
One group of papers [§, H 0, §J H from the very beginning 
describes the system as consisting of two (interacting) frac- 
tions: the condensate and the thermal cloud. This ambitious 
program is very demanding numerically. The other group in- 
terprets the high energy solutions of the time-dependent GP 
equation as describing the full condensate plus thermal cloud 
system [ 10 1. As a rule, these authors were able to identify the 
condensate as a definite part of the system only in the some- 
what academic case of the gas in a rectangular box with peri- 
odic boundary conditions. In this case the condensate is just 
the zero momentum component of the wave function. The aim 
of this Rapid Communication is to show that much more may 
be achieved with the readily available high-energy solutions of 
the GP equation in harmonic traps. It is also possible to split 
this solution into a sum of the condensed and uncondensed 
parts, define the condensate wave function and approximately 
estimate the temperature of the resulting system. 

The Hamiltonian H of the system takes the form: 
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[4*(r, f), (r', f)] = <5(r - r')- The first term describes par- 
ticles with mass m trapped in a potential of a spherically 
symmetric harmonic oscillator of frequency lo while the sec- 
ond term describes two-body interactions. Here we have as- 
sumed that particles interact via a contact potential V(r - r') = 
AnTi 2 a s 5(r - r')/m, where a s is the s-wave scattering length. 
The Heisenberg equation originating from this Hamiltonian 
acquires the following form 
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A full operator solution of the non-linear Eq.(Q) is not known. 
In addition, the identification of a condensate phase is also a 
subtle issue. The only exception is a system of particles in a 
periodic box. Here the symmetry of the problem helps: nat- 
ural eigenmodes of the system are plane waves with a quan- 
tized momentum k = 2n( j\ , ji, jj,)/L (y,- being integer, L - box 
size). Therefore, the field operator can be expanded in these 
modes: 



(3) 



where *P(r, f) is a field operator that destroys a particle at po- 
sition r and obeys standard bosonic commutation relations 



A Bose-Einstein condensate can be uniquely associated with 
the zero momentum mode, k = 0. The annihilation and cre- 
ation operators a^, a R satisfy non-linear equations following 
fromEq.(|). 

Different approximate approaches mentioned above [|5j |[ 
0, 0, |J] have been introduced in order to solve the problem. 
Here we want to utilize a method which has been extensively 
and successfully explored in quantum optics. A quantum field 
operator describing a coherent electromagnetic field (such as 
laser light) can be replaced by a complex valued classical 
field. Such a substitution is justified for these modes that are 
highly occupied. Only then are quantum fluctuations negligi- 
ble and the non-vanishing commutator may be ignored. This 
kind of approximation has been used recently in JlO[ ] for the 
matter field: — > a^. Note that £ k l fl kl 2 = N and the con- 
densate occupation is equal to no = l«ol 2 - On the other hand, 
it follows from Eqs.(||)-(Q) that the semi-classical approxima- 
tion is equivalent to a substitution: *T(r, t) — > VAf*T(r, t) - 
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-jjii 2k a k exp(-/k • r), where *F(r, t) fulfills the standard GP 
equation. According to conventional wisdom, this equation 
describes a pure condensate at zero temperature, hence the 
whole classical field *P(r, t) represents the condensate popu- 
lated by n = N particles [|nj 




FIG. 1: Cross-section of the instantaneous snapshot of the density 
distribution in the z = plane of the Bose gas with the total energy 
E = 70 %u). In this and all the following figures lengths are given in 
units of d = \/Ti/mti>. 

It seems that the two contradicting interpretations of the 
high-energy solution coexist in the literature. In the follow- 
ing we will clarify this seeming incongruity. We focus our 
attention on experimentally relevant system of 100 000 87 Rb 
atoms trapped in the spherically symmetric harmonic poten- 
tial of frequency u> = 27rl00Hz. Our procedure is as follows. 
First, we define the wave function of the system on a spatial 
three-dimensional grid. Initial values of the wave function at 
each point are chosen in accordance with the constraints of 
a fixed energy and particle number. Next, we propagate this 
state in time, solving the time dependent GP equation with 
the help of the Fast Fourier Transform. It occurs that for the 
given total energy and particle number the same steady state is 
reached after several milliseconds. We have checked that the 
steady state attained does not depend on the choice of the ini- 
tial wave function (in particular, on the fact whether its phase 
is random or not) as long as the number of particles and the 
total energy are kept constant. 

The high-energy solutions of the GP equation have strik- 
ing features [ |ll|. A snapshot of the density distribution 
is shown in Fig.|l[ The density is extremely irregular and 
exhibits a number of sharp spikes changing its shapes and 
positions on very short time scales. Typical methods used 
to monitor trapped atomic condensates involve optical tech- 
niques. A condensate is illuminated by laser light which, af- 
ter passing through the atomic system (and some optical ele- 



ments), is monitored by a CCD camera. The exposure time At 
varies from a few microseconds to hundreds of milliseconds 
and within this time the condensate density undergoes rapid 
changes due to the fast dynamics but also because of quantum 
fluctuations. The nature of the observation process introduces 
a kind of smearing. In fact, in such a long exposure a time- 
averaged rather then an instantaneous single-particle density 
is monitored. The existence of short time scales in non-linear 
many-body dynamics has to be taken into account. The object 
of physical significance is therefore a time-averaged single- 
particle density matrix: 
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This coarse graining procedure destroys the purity of a state 
of the system. After the time averaging all irregularities of the 
density are being smoothed out especially in the inner part of 
the density profile (see Fig^J). The resemblance to the well 
publicized photographs of the experimental condensates at in- 
termediate temperatures is striking jl3[]. Varying the energy 
of the gas we can scan the whole range from the pure conden- 
sate all the way to a nearly critical condition of no condensate. 
The sample profiles of the time-averaged column density (in- 
tegrated along the z-axis) are shown in Fig^. In parts A and 
B of Fig.|| note the bimodal structure of the distribution: the 
central peak corresponding to the Bose-Einstein condensate 
and the broad background identified with the thermal cloud. 

We can provide a quantitative analysis of the resulting av- 
eraged state of the Bose gas. To this end we recall a classical 
definition of the condensate by means of the spectral decom- 
position of the single -particle density matrix [|l4]]. Diagonal - 
ization of a single-particle time-averaged density matrix leads 
to natural eigenmodes t/r,(r, f) and eigenvalues rij/N: 

p(ri,r 2 ,f) = J] (ri,f)fcta, 0- (5) 
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The system may be viewed as a mixture of many coupled co- 
herent modes i^, whose occupation is n,. The self-consistency 
of the model requires that population of each single mode 
is large, «, > 1. Otherwise the substitution of the field op- 
erator T^r, t) by a wave function l F(r, t) would not be justi- 
fied. The existence of the dominant eigenvalue hq/N of the 
order of unity signifies the presence of the condensate with 
the corresponding eigenvector iffg being its wave function. As 
the modes and their occupation are known only after the time 
averaging, the verification of the semi-classical criterion can 
only be done a posteriori. We try a number of grids and 
choose the one that yields the highest occupation of all modes, 
justifying the semi-classical approximation. 

The time averaging followed by a diagonalization of the 
density matrix gives the spatial modes of the system. These 
modes are typically not known even in a stationary case. Still, 
the only exception is the box with periodic boundary condi- 
tions. Let us illustrate our method using again this simple 
case. 



FIG. 2: Time-averaged stationary column density distribution for 
a Bose gas plotted for different values of total energy E. A E = 
70.00 ha> (30% atoms in condensate). B E = 23.67 hoi (72% of atoms 
in condensate). C E = 13.04S<y (pure condensate). 



Substitution of all annihilation operators by complex ampli- 
tudes 2k — > at gives the following expression for the single- 
particle density matrix: 

p(ri,r 2 ,0 ~ - ^ exp[/k • (n - r 2 )]|a k (f)| 2 (6) 

k 

and the time averaging leaves only diagonal elements k = k' 
as the off-diagonal ones oscillate rapidly and get rapidly de- 
phased, a(t)^a(t) k > ~ <5k,k'l a k(OI 2 - Obviously, the |a k | 2 are pop- 
ulations of different modes. Only if one introduces a kind of 
coarse graining leading to a suppression of the off-diagonal 
elements of the single-particle density matrix can one identify 
the k = mode as a Bose-Einstein condensate with an occu- 
pation given by |aol 2 ■ Without this additional assumption the 
whole complex field *P(r, f) describes one coherent, dynami- 
cally evolving Bose-Einstein condensate without any thermal 
cloud. This kind of interpretation is used in JTl| , |T2[ ]. On the 
contrary, identification of the k = momentum component of 
^(r, f) with a Bose condensate has been directly assumed in 
JTc|] . Our analysis solves this apparent contradiction. By ex- 
amining a detection process and the relevant time scale we can 
uniquely determine the condensate fraction, its wave function 
as well as the structure of excited modes of the interacting 
system. Moreover, the method is no longer limited to an aca- 
demic problem of the uniform system. A steady-state single- 
particle density matrix for the spherically symmetric trapping 
potential must have eigenvectors that are simultaneously di- 
agonalizing the angular momentum operators. Therefore, the 
eigenvectors must be proportional to the spherical harmonics 
Y[ m . What remains is the one dimensional diagonalization in 
the radial variable of the following projection of the density 
matrix: 



Plm(r, r') = J p(r, Q; r', Of) F* (Q) F, m (Q') dCl dtl' , (7) 

where the integration is performed over solid angles Q. and Q.' 
associated with the corresponding particle coordinates. We 



expect the condensate to be present in the zero angular mo- 
mentum component of the single-particle density matrix. In- 
deed, in the inset of Fig.|| we show a typical distribution of the 
eigenvalues of the / = part of the density matrix with one 
dominant eigenvalue and with the corresponding eigenfunc- 
tion plotted in Fig.|| together with the whole time-averaged 
radial density distribution. 
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FIG. 3: Time-averaged radial density profile of a stationary state of 
a Bose gas. The total energy is E = 70.00 hu>. Density of the whole 
system - solid line; profile of the Bose-Einstein condensate (defined 
as a dominant eigenmode of poo(n> r i)) - dashed line. All eigenval- 
ues of poo(n, ri) are shown in the inset. 

The last point raised here is the question of temperature. In 
principle we could do something similar to a typical experi- 
ment: switch off the trapping potential, let the gas expand and 
analyze the properties of the thermal part. Numerical limita- 
tions and the constraining condition of the high occupation of 
each eigenmode make it hard. Alternatively, we can follow 
the procedure of [|l5|| and fit the profile of the thermal fraction 
of the ideal Bose gas to the outer part of our averaged density. 
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This yields a reasonable, order of magnitude result. The tem- 
perature dependence of the condensed fraction together with 
its large uncertainty is shown in Fig^. In the inset we plotted 
the energy dependence of the number of atoms in the con- 
densate. We compare our results with the ideal gas and with 
a two-gas estimate (finite-temperature Hartree-Fock scheme 
fll6|]). We see a reasonably good agreement between all three 
calculations. 
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FIG. 4: Condensate fraction versus temperature or (the inset) the to- 
tal energy of the system. The solid line represents the ideal gas in the 
thermodynamic limit, the dashed line depicts the results of jl6| ob- 
tained within a finite-temperature Hartree-Fock scheme, while dots 
show our data. Error bars indicate a range of temperatures for which 
we obtain acceptable fits of the ideal Bose gas distribution to the 
outer wings of the thermal cloud. T c is the ideal gas critical temper- 
ature in the thermodynamic limit and E = ha>. 

To summarize: we have shown a simple way of numeri- 
cally simulating the stationary dynamics of a weakly inter- 
acting trapped Bose gas at finite temperature, which uses a 
semi-classical representation of the matter field and, in ac- 
cordance with the experimental procedures, stresses the role 
of a finite exposure time when photographing the condensate. 
In the self-consistent determination of the condensate fraction 
we rely solely on the classic criterion of Onsager and Penrose. 
This way we avoid an arbitrary splitting of the system into the 
condensed and thermal components from the very beginning. 
The method may be used to model the impact of thermal fluc- 
tuations on the dynamical processes with the condensate, such 
as solitons or vortices. The next step in going beyond the ap- 
proximations employed in our model is the estimation of the 
influence of quantum corrections to the classical fields, in par- 
ticular the study of the corresponding time scales. 
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